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Nonlinear oscillators can mutually synchronize when they are driven by common external impulses. 
Two important scenarios are (i) synchronization resulting from phase locking of each oscillator to 
regular periodic impulses and (ii) noise-induced synchronization caused by Poisson random impulses, 
but their difference has not been fully quantified. Here we analyze a pair of uncoupled oscillators 
subject to common random impulses with gamma-distributed intervals, which can be smoothly 
interpolated between regular periodic and random Poisson impulses. Their dynamics are charac- 
terized by phase distributions, frequency detuning, Lyapunov exponents, and information-theoretic 
measures, which clearly reveal the differences between the two synchronization scenarios. 



I. INTRODUCTION 



Rhythmic elements are basic building blocks of Nature at human scales. They play particularly important roles in 
the functioning of biological systems, such as cardiac muscle cells, pacemaker neurons, and animals or plants following 
circadian or circannual rhythms. Recently, synchronization among non-interacting rhythmic elements induced by 
common fluctuating external force has attracted much attention in diverse fields. It has been demonstrated for a wide 
class of rhythmic elements ranging from lasers and electronic circuits to biological systems [U-IH . For example, reliable 
synchronous response of neurons due to common or shared input signals has been actively discussed in neuroscience . 
In ecology, it is known that, due to common climate fluctuations, populations of plants exhibit large-scale synchronized 
flowering and production of seed crops also fluctuate synchronously from year to year [Iol - [T2j ]. and such phenomena 
are generally termed the "Moran effect" [l3| . 

We often characterize a rhythmic element as a limit-cycle oscillator. Populations of coupled limit-cycle oscillators 
show a variety of interesting collective behavior, including mutual synchronization, wave propagation, and chaos (l4l - 
It is also well known that uncoupled limit-cycle oscillators can mutually synchronize when they are driven by 
common impulses pit IbH. [l6l - [l9j . The simplest classical situation is, of course, synchronization due to phase locking 
of an oscillator to common periodic impulses [l7| . Another interesting situation is noise-induced synchronization as 
mentioned above, caused e.g. by common random Poisson impulses |3l.ll6l.ll8l.ll9||. The oscillators can also synchronize 
when the driving impulses have intermediate statistics between the periodic and the Poisson impulses, as shown by 
Yamanobe et al. [l8| for a model of pacemaker neurons driven by gamma-distributed impulses. 

This prompts the question, what is the difference between the synchronization due to phase locking and noise- 
induced synchronization? Though drive-response configuration of the impulse source and the oscillators is the same 
for both cases, it is expected that there should be some difference in their synchronization dynamics, reflecting 
different characteristics of the driving impulses. In this paper, we address this issue by considering uncoupled limit- 
cycle oscillators driven by gamma impulses, which can be smoothly interpolated between regular periodic and random 
Poisson impulses. We examine the transition from phase locking to noise-induced synchronization as the statistics 
of the gamma impulses is varied, and quantify their difference using phase distributions, Lyapunov exponents, and 
information-theoretic measures . 

Effect of common gamma impulses on limit-cycle oscillators was previously treated in the paper by Yamanobe 
et al. [la ], but its main focus was not on the transition between the two types of synchronization but rather on 
physiologically realistic characterization of pacemaker neurons. We here conduct a systematic, quantitative analysis of 
the transition between the two different synchronization dynamics. We will reveal that the two types of synchronization 
dynamics are clearly different in many aspects, e.g. their stability, fluctuations, and statistical dependence on the 
driving impulses. 

This paper is organized as follows. In Sec. II, we introduce a model of uncoupled oscillators subject to common 
gamma impulses and demonstrate its synchronization dynamics for different types of driving impulses. In Sec. Ill, we 
analyze phase distributions of the oscillators using Frobcnius-Perron-type equations. In Sec. IV, we focus on frequency 
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detuning of the oscillators due to impulsive driving. In Sec. V, we examine Lyapunov exponents and their fluctuations 
in the synchronized states. In Sec. VI, we characterize mutual relationships among the impulse source and the driven 
oscillators by utilizing information-theoretic measures. Section VII summarizes the results. 

II. LIMIT-CYCLE OSCILLATORS DRIVEN BY GAMMA IMPULSES 

In this section, we introduce phase oscillators driven by common gamma impulses and demonstrate their synchro- 
nization dynamics. 

A. Random phase map 

Based on the standard phase reduction theory of limit cycles pH [I3 - [l6l Il9j. we can describe a limit-cycle oscillator 
using only its phase variable, 9 € [0, 1), defined along its stable limit cycle provided that the intervals between the 
impulses are sufficiently longer than the amplitude relaxation time of the oscillator state to the limit-cycle orbit. The 
key quantity for this description is the phase response curve (PRC) G{9) of the oscillator [3], which encapsulates 
its dynamical properties. G{9) is a periodic function on [0, 1) that gives asymptotic response of the oscillator phase 
to an impulsive perturbation given at phase 6. It has been measured experimentally, e.g., for repetitively firing 
neurons [2(| [2l[. When the limit-cycle oscillator is kicked by an impulse at phase 6, the phase is mapped to a new 
value F(6) = 9 + G(9). This function F(9) is sometimes called the "phase transition curve" in the literature. 

The dynamics of the oscillator between two successive impulses consists of a jump caused by the first impulse and 
subsequent free rotation until the arrival of the second impulse 0, [lij [22[ ■ The oscillator is also subject to small 
temporal fluctuations of various origins. We assume that each oscillator receives common driving impulses at times 
{ti, t2, ■ ■ ■ , t n , ■ ■ ■ }, and denote the phase of the oscillator just before the nth impulse by 9 n £ [0, 1). The effect of the 
nth impulse is to map the oscillator phase 9 n to F(9 n ). To incorporate the effect of small fluctuations, we also apply 
a weak additive independent noise rj n before the mapping. The phase 4> n just after the nth impulse is thus given by 
<t>„ = F(9 n + Tin). After receiving the nth impulse, the oscillator rotates with a constant frequency until t n +i at which 
the (n + l)th impulse arrives. Therefore, the phase 9 n +i just before the (n + l)th impulse is given by 

9 n+1 = F(9 n + r) n ) + ujt 71 mod 1, (1) 

where r„ = t n+ i — t n is the inter-impulse interval, u is the frequency of the oscillator, and we have taken modulo 1 of 
the phase to confine it in [0, 1). Since r„ and rj n are random variables, this equation gives a random map, which we will 
call a "random phase map" hereafter. When we consider multiple oscillators, r„ is common to all oscillators, whereas 
rjn is different from oscillator to oscillator. Note that r) n represents small temporal fluctuations of the dynamics of 
oscillators, but not static heterogeneities in their natural frequencies. Mean frequency of the driven oscillator exhibits 
qualitatively different dependence on the driving impulse between phase locking and noise-induced synchronization, 
as we show in Sec. IV. 

As the simplest and typical example, we assume that the PRC is given by a sinusoidal function G(9) = csin(27r#) 
in the following numerical analysis, so that the phase map is given by 

F(6) =0 + csin(27T0). (2) 

The parameter c controls the magnitude of the nonlinearity, which may also be regarded as the intensity of the 
impulse. Such a PRC with both positive and negative lobes are called Type-II (8l. [23l - [25j . which generally arises near 
the Hopf bifurcation of limit-cycle oscillators [J I23W27I ] . Another typical example is the Type-I PRC with a single 
positive lobe, e.g. G(9) = cjl — cos(27r6>)|, which generally arises near the SNIC (saddle-node on invariant circle) 
bifurcation of oscillators [1, [25|, [28|. However, difference between these two PRCs can roughly be eliminated by 
simply shifting the frequency lo and the origin of the phase 9 in the present setup, and qualitatively similar results are 
expected for both PRCs. Actually, the two types of PRCs yielded almost the same results in our numerical analysis. 
Thus, the type-II PRC, Eq. ©, gives us general insights and we only present the results for this case hereafter. 

B. Gamma impulses 

We consider driving impulses that have "intermediate" statistics between regular periodic impulses and random 
Poisson impulses. As such impulses, we adopt the gamma impulses [l8j . whose inter-impulse interval r obeys the 
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FIG. 1: (Color online) The gamma distribution W(t; a, b) for (a) a = 1, b ■■ 
The mean value is fixed at (t) = ab = 1. 



1, (b) a = 30, b = 1/30, and (c) a = 1000, fe = 1/1000. 



gamma distribution [T8|, [29| , 

where a and 6 are real positive parameters. The mean interval is given by (r) = r o °° W(t; a, &)rdr = ab. The gamma 
distribution has the following properties. Firstly, when a = 1, W-^r; a, b) gives an exponential distribution, 

W{r;a = l,b) = ^-, (4) 

which means that the impulses obey the standard Poisson random process. Secondly, by taking the limit a — > oo and 
b — > with (r) = ab fixed, the gamma distribution approaches the delta function, 

W(r;a,b) a ^ 5(r-(r)). (5) 

In this limit, the inter-impulse interval is always (r), so that the impulses become periodic. 

Thus, the parameter a, which we call a shape parameter, determines the property of the gamma impulses as shown 
in Fig.[TJ By varying the value of a between 1 and oo with (r) = ab fixed, the gamma impulses can exhibit intermediate 
properties between random Poisson and regular periodic impulses while keeping the same mean inter-impulse interval. 
We examine uncoupled oscillators driven by gamma impulses and see how synchronization dynamics of the oscillators 
depends on the shape parameter a in the following. 

The gamma distribution has several nice mathematical properties and naturally arises under a few simple assump- 
tions when the mean and the irregularity of impulse sequences, which respectively correspond to the mean interval (r) 
and the shape parameter a [29j ]. are given. We thus use the gamma impulses in the present paper, though alternative 
approaches for varying the regularity of the driving signal should also be possible, e.g., by using chaotic dynamical 
systems [30| . 



C. Synchronization by common driving impulses 

We here demonstrate synchronization of uncoupled oscillators induced by common impulses for the periodic, Poisson, 
and intermediate cases by direct numerical simulations. We fix the magnitude of nonlinearity (or the intensity of 
impulses) of Eq. ((2J at c = 0.1 throughout our numerical simulations, which is sufficiently small and therefore the 
map is not chaotic even for strictly periodic impulses. This is because we are focusing on the synchronization dynamics 
of limit cycles (note that the sinusoidal Type II phase map with periodic impulses is nothing but the well-known circle 
map fl6|). Strong Poisson impulses may also yield positive Lyapunov exponents and lead to desynchronization of 
oscillators as shown in Ref. 0, but we do not consider such a situation in the present paper. 

Figure [2] shows typical synchronization dynamics of uncoupled oscillators induced by three types of common im- 
pulses, where zero-crossing points of 10 oscillators are shown in raster plots. The mean interval of the impulses is set 
at (t) = ab = 1. The period of the oscillator is also taken asT = l/oj = l and is thus equals to the mean interval. 

(i) Phase locking to periodic impulses [Fig. [2] (a)]. If the impulses are periodic and their intervals are nearly equal to 
(or more generally rational multiples of) the intrinsic rotation period of the oscillators, namely, if they are resonant, 
each oscillator becomes phase locked to the impulses. As a consequence, uncoupled oscillators driven by common 
periodic impulses synchronize with each other. 
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(ii) Noise- induced synchronization by Poisson impulses [Fig. [2](b)]. The oscillators also synchronize when they are 
driven by common Poisson impulses of appropriate intensity. It has been theoretically and experimentally shown that 
uncoupled limit-cycle oscillators subj ect to weak common Poisson impulses generally synchronize with each other, 
irrespective of their details 0, EE EH- 

(hi) Synchronization induced by gamma impulses [Fig. [2(c)]. Uncoupled oscillators subject to the gamma impulses 
with intermediate statistics between the periodic and Poisson impulses can also synchronize mutually (l8| . 

Thus, the uncoupled oscillators synchronized with each other by the effect of the common impulses for all values of 
the shape parameter. It is however not easy to catch the quantitative differences in their synchronization dynamics 
just from these figures. In the following sections, we will characterize the differences for varying types of common 
impulses using phase distributions, frequency detuning, Lyapunov exponents, and information-theoretic measures. 
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FIG. 2: (Color online) Synchronization of uncoupled oscillators induced by common impulses. Rasters show zero-crossing events 
of individual oscillators. Each arrow shows the time when an impulse arrives, (a) Phase locking (a — 1000). (b) Noise-induced 
synchronization (a = 1). (c) Synchronization induced by gamma impulses (a — 30). 



III. PHASE DISTRIBUTIONS 

In this section, we introduce Frobenius-Perron-type equations 0, |U H2| that describe evolution of phase distribu- 
tions. Using them, we examine how the relations among the oscillator phases and the inter-impulse intervals depend 
on the shape parameter, which reveal differences between phase locking and noise-induced synchronization. 
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A. Single-oscillator Frobenius-Perron equation 



Let us consider a single-oscillator probability density function (PDF) P n {9) of the phase 6 n at time step n, just before 
the nth impulse arrives. To obtain a Frobenius-Perron equation describing the dynamics of P n (6), it is convenient to 
consider a PDF P n (9) of the "unwrapped" phase 9„ defined in (—00,00), which obeys the following random phase 
map without the modulo 1, 

0„+i = F(6 n +r) n )+LOT n . (6) 
The PDF P n+ i(9) of 9 n+ i is given by the following Frobenius-Perron equation (32l. |33|: 

/>1 />oo />oo 

P n +i(S)= d9 dr dr)P n (6)W(T)R(ri)S(6-F(e + r))-UT), (7) 

JO JO J-00 

where W(r) is the PDF of the inter-impulse intervals, namely, the gamma distribution given in Eq. ([3]), and 

R(n) = -L= cxp (-f-) (8) 



V2^2 V 2a 2 

is the probability density function of the independent additive noise, which we assume to be zero-mean Gaussian of 
variance a 2 . The probability density function of the wrapped phase, 9 = 8 mod 1 € [0, 1), can be calculated as [34j 

00 

Pn{0)= P n(Q + k). (9) 

k— — 00 

Thus, from Eq. (O, the Frobenius-Perron equation describing P n (9) is obtained as 

00 pi poo poo 

P n+ i{6)= V dO dr dnP n (e)W(T)R(r))5(6-F(8 + ri)-u>T-k), (10) 

k=-oo J ° J ° J -°° 

which describes the single-oscillator PDF P n {9) of the phase 9 n just before the nth impulse. 

Stationary solutions of Eq. (|T0|) gives the PDF of the oscillator phase driven by gamma impulses and additive 
noise in statistically steady states. In Ref. [l9j], we have treated a similar Frobenius-Perron equation perturbatively 
to obtain the stationary phase PDF of the oscillator under weak Poisson impulses and calculated the Lyapunov 
exponents. In Ref. (3lj . Doi et al. analyzed the spectral properties of a similar Frobenius-Perron equation (driven by 
periodic impulses) to characterize noisy phase locking of a limit-cycle oscillator. 

B. Two-oscillator Frobenius-Perron equation 

To investigate phase relationships between a pair of oscillators (denoted here as A and B) subject to common 
impulses, we introduce a joint PDF P n {9 A , 9 B ) of their phases 9 A and 9 B just before the nth impulse. The dynamics 
of the pair is given by 

<+l = F(0n + Vn) + UT n mod 1, 

Cl = F (°n + Vn) + "T n mod 1, (11) 

where r\ A and are independent additive noise terms (r„ is common to both oscillators). Similarly to the single- 
oscillator case, we can derive a Frobenius-Perron equation describing the two-oscillator phase PDF P n (9 A ,9 B ) as 

P n+1 (9 A ,9 B )= J2 E d9 A ' d9 B ' dr dn A ' drj B ' P n (9 A ' ,9 B ' )W(t)x 

R(n A ')R( V B ')S (V - F(9 A ' + r) A ') -lot- k A ) 5 [e B - F(9 B ' + V B ') - lot - k B ) ■ (12) 

The stationary solution of Eq. p^|) gives a PDF of the pair of oscillator phases driven by common impulses and 
independent additive noise in the statistically steady state. 

In Ref. [22| , we have analyzed a similar two-oscillator Frobenius-Perron equation by averaging out the fast dynamics 
of the mean phase to obtain an approximate Frobenius-Perron equation describing only the phase difference of two 
oscillators driven by common Poisson impulses. Since we arc interested in pair phase distributions, not only in the 
distribution of phase differences, we directly solve Eq. (|12p numerically for gamma impulses in the following. 




FIG. 3: (Color online) Stationary single-oscillator probability density function P(9) of the phase 9 just before the arrival of 
each impulse. Blue curves are obtained by numerically solving the Frobenius-Perron equation and red squares are the results 
of direct numerical simulations of the random phase map. (a) Noise- induced synchronization (a = 1). (b) Synchronization 
induced by gamma impulses (a = 30). (c) Phase locking (a = 1000). 
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FIG. 4: (Color online) Stationary two-oscillator joint probability density function P(6 ,8 B ) of the phase pair 9 and 9 B just 



before the impulses. I 
(c) Phase locking (a 
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1000). 



a — 1). (b) Synchronization induced by intermediate impulses (a = 30). 



C. Numerical results 



We here present stationary phase PDFs obtained by numerically solving the Frobcnius-Pcrron equations. We set 
the frequency of each oscillator at lu — 1. The period of the oscillator is T = 1/lo = 1 and is equal to the mean 
period of impulses (r) = 1. Thus, when the shape parameter a is sufficiently large and the impulses are nearly 
periodic, phase locking should occur. In the numerical analysis, each variable is discretized using 100-200 grid points. 
The calculation cost can be drastically reduced by devising numerical algorithms that use Fourier representation in 
calculating convolutions of the Frobenius-Perron equations (see Appendix for details) . 




FIG. 5: (Color online) Stationary joint probability distribution P(r, 9) of the inter-impulse interval r and the phase 9 imme- 
diately after this interval, (a) Noise-induced synchronization (a — 1). (b) Synchronization induced by intermediate impulses 
(a = 30). (c) Phase locking (a = 1000). 
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1. Single-oscillator phase distributions 

Figures [3ja)-(c) show the stationary PDF P(9) of a single oscillator phase 9 just before each impulse for nearly 
periodic (a = 1000), intermediate (a = 30), and Poisson (a = 1) cases. Numerical solutions of the Frobenius-Perron 
equation (|10[) and the results of direct numerical simulations of the random phase map (TTJ) agree nicely. The intensity 
of the weak independent Gaussian noise is a = 0.01. In all cases, the oscillators synchronize with each other (to the 
independent noise level). However, depending on the shape parameter a, the phase PDF P(9) significantly differs. 
For nearly periodic impulses, P(9) has a sharp peak, indicating that the oscillator phase is actually entrained to the 
driving impulse with some fixed phase shift. In contrast, P(8) is roughly uniformly distributed for Poisson impulses, 
which means that the relationship between the oscillator phase and the impulse timing is not fixed and hence not 
entrained. In the intermediate case, P(0) has a broad but still clear one-humped shape, implying that the phase is 
not completely entrained to the impulses but still possesses a certain level of coherence with respect to the driving 
impulses. 

2. Two-oscillator phase distributions 

Figures Hta)-(c) show the stationary joint phase PDF P(9 A ,9 B ) of a pair of oscillators for nearly periodic (a = 
1000), intermediate (a = 30), and Poisson (a = 1) impulses obtained by numerically solving the Frobenius-Perron 
equation (|12p. In all cases, the oscillators synchronize with each other, so that their phases are distributed along the 
diagonal line, 9 A = 9 B . However, their distribution strongly depends on the shape parameter a. For nearly periodic 
impulses (a = 1000), the oscillator is almost phase- locked to the impulses and thus the PDF has a sharp peak near 
the center, indicating that the oscillators are not only mutually synchronized but both of them are entrained to the 
impulses. As the parameter a decreases, the distribution becomes broader, and for Poisson impulses (a = 1), the 
phases 9 A and 9 B are broadly distributed along the diagonal line, indicating that they are synchronized but not 
entrained by the impulses anymore. 

3. Joint distributions of the inter-impulse intervals and the oscillator phases 

To analyze how the driving impulses affect the oscillator phase, we also calculate the joint PDF of the inter- 
impulse interval and the oscillator phase just after this interval in the statistically steady states. The stationary 
joint PDFs P(r, 9) = P(9\t)W(t) of the impulse interval r and the phase 9 observed immediately after this interval 
are calculated from the Frobenius-Perron equation. Figures [5ja)-(c) show stationary joint PDFs P(t,9) for nearly 
periodic (a = 1000), intermediate (a = 30), and Poisson (a = 1) cases. The distribution has a sharp peak in the 
nearly phase-locked case (a — 1000), indicating that the inter- impulse interval and the oscillator phase just after this 
interval have almost a one-to-one correspondence, namely, the oscillator phase is almost entrained by the impulses. 
As we decrease the shape parameter a, the distribution gradually elongates, and in the Poisson limit (a = 1), such 
clear correspondence is lost (but they still retain a certain degree of correlation). 

Thus, all phase distributions clearly capture the essential difference between the phase locking and the noise- 
induced synchronization. The relation among the oscillator phases and the inter-impulse intervals significantly differs 
depending on the shape parameter a, even though the oscillators themselves are always mutually synchronized. As a 
is increased, the oscillators tend to be more strictly phase locked to the driving impulses, whereas for smaller a, their 
dependence becomes weaker. This observation will be quantified by information-theoretic measures in Sec. VI. 

IV. FREQUENCY DETUNING 

In this section, we consider how the dynamics of the oscillator is modulated by the driving impulses. Specifically, we 
examine the frequency detuning i.e., the difference between the mean frequency of the driven oscillator and that 
of the driving impulses, for varying values of the shape parameter. This analysis will reveal another clear difference 
between the two types of synchronization scenarios. 
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A. Mean frequency of the driven oscillator 

From the phase PDF P(9) and the PRC G(9), we can estimate the mean frequency of the impulse-driven oscillator 
as follows. Using the random phase map for unwrapped phase 9 n , Eq. ©, the oscillator phase just before the iVth 
impulse is given by 



N-l 

?jv = 0i + 



J2 (G(0i + m) + wr<) . (13) 



Therefore, long-time mean frequency u/ of the driven oscillator can be calculated from this equation as 



uj = lim 



jV->oo N — 1 

/ 1 N-l 1 N-l ' 

= lim G(9 l + r>i) H 

\ i=l i=l 

/>1 />00 

= / d6P{9)G(9) +u drW(r) 



d9P(9)G(9) +u, (14) 

where we have replaced the long-time average by a statistical average and used the periodicity of the PRC G(9) and 
the normalization condition J °° dTW(r) = 1. Thus, the mean frequency detuning between the oscillator and the 
impulses is given by 

u'-Cl= [ d9P(9)G(9) + uj -n, (15) 
Jo 

where f2 = l/(r) is the mean frequency of the driving impulses. 



B. Numerical results 



We fix the mean frequency of the driving impulses at f2 = l/(r) = 1 and vary the shape parameter a and the 
natural frequency w of the oscillator in the numerical simulations. Figure |6ja) plots the frequency detuning u)' — Q 
on the (a, u>) plane. For large a and nearly resonant natural frequency, u ~ fl = 1, the conditions for phase locking 
are satisfied and therefore the mean frequency uj' of the oscillator is locked to that of the driving impulses f2, yielding 
a plateau in the figure. In contrast, for small a, the mean frequency u>' never locks to f2. 

Figure EJb) plots u>' — fl as a function of u> at a = 1000, i.e., for approximately periodic impulses. Results of direct 
numerical simulations agree well with the theoretical estimate, Eq. ([TS")) . with the PDF P(9) obtained numerically 
from the Frobenius-Perron equation. A plateau satisfying uj' = can clearly be seen near uj = f2. In the intermediate 
case (a = 30) shown in Fig. Etc), clear phase- locking region no longer exist, however, the frequency of the oscillator 
is modulated by the driving impulses. In the Poisson case (a — 1), the mean frequency of the driven oscillator is not 
affected by the driving impulse as shown in Fig. Htd). 

This observation provides us with further evidence on the difference between the two synchronization scenarios. 



Actually, as pointed out by Yoshimura et al. in Rcf. [351 ] . absence of frequency locking between the oscillators is 



characteristic to common-noise-induced synchronization, in sharp contrast to ordinary phase locking due to mutual 
coupling (though they considered weak Gaussian noise instead of random impulses in their work) . 



V. STABILITY OF SYNCHRONIZED STATES 



In this section, to characterize the synchronization dynamics, we focus on Lyapunov exponents and their fluctuations 
in the synchronized states. As we will sec, differences between synchronization scenarios are well characterized by the 
fluctuations of the Lyapunov exponent. 
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FIG. 6: (Color online) Frequency detuning, (a) Dependence of the frequency detuning uJ — on the shape parameter a and on 
the natural frequency ui of the driven oscillator, measured by direct numerical simulations of the random phase map. (b)-(d) 
Dependence of the frequency detuning on the natural frequency cu of oscillators for several values of the shape parameter a. 
Blue solid lines are the results of the integral calculation, Eq. (|15[) . and red squares are obtained by direct numerical simulations 
of the random phase map. (b) Phase locking (a = 1000). (c) Synchronization induced by gamma impulses (a = 30). (d) 
Noise- induced synchronization (a = 1). 



A. Lyapunov exponents and its variance 

To quantify statistical linear stability of the synchronized states, we calculate the Lyapunov exponent and its 
variance. Let us consider a pair of oscillators and denote their phases at time n as 9 n and 0' n , respectively. Linearized 
evolution of a small phase difference A„ = 6' n — 6 n is given by 



A n+1 =F'(8 n + T ln )A r , 



(16) 



where F'(6) = dF(6)/dO is the instantaneous linear growth rate of the phase difference. The deviation Ajv at large 
time step N is thus given by 



Ao 



JV-l 



I |F'(0„+77„)|^exp[iV(A)], 



(17) 



71=0 



where we have introduced the Lyapunov exponent (A) , defined by a long-time average of the linear growth rates as 



(A) = lim — In 



N-l 



d6 



dr)P(6)R(r]) In \F' (6 + rf)\ . 



(18) 



In the last expression, we have replaced the long-time average by a statistical average over the stationary phase PDF 
P{ff) and the PDF R(rj) of the independent noise. The synchronized state is stable if (A) is negative. 
Similarly, variance of the linear growth rates, var(A) = (A 2 ) — (A) 2 , can be calculated as 



var(A) = / d6 I dj]P(e)R(rj){\n\F' (9 + } 2 

— oo 

d6 f dr)P(e)R(r))ln\F' (6 + ri)\\ . (19) 

As we will see, even if the Lyapunov exponent (A) takes the same value, its variance var(A) can significantly differ 
depending on the shape parameter a. 
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FIG. 7: (Color online) Dependence of the Lyapunov exponent A on the shape parameter a and on the frequency w of oscillators, 
(a) Mean Lyapunov exponent (A) . (b) Variance of Lyapunov exponent var(A) . (c) Ratio of the mean to the variance, (A) /var( A) . 
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FIG. 8: (Color online) Dependence of the mean and variance of Lyapunov exponent A on the parameter a. The curves are 
calculated from the Frobenius-Perron equation and the symbols are the results of direct numerical simulations of the random 
phase map. (a) Mean (A), (b) Variance var(A). 



B. Numerical results 



Figure [TJa) plots the Lyapunov exponent (A) on the (a,ui) plane, showing its dependence on the shape parameter 
a and on the natural frequency of the oscillator u>. The mean period of the impulses is fixed at (r) = 1. For 
sufficiently large a, there exists bell-shaped regions near the resonant frequencies (i.e., u — 1 and u = 2), where 
the Lyapunov exponent takes large negative values. These regions correspond to nearly phase-locked states and 
may be seen as a kind of Arnold tongues. The Lyapunov exponent also takes small negative values outside this 
region, which corresponds to the noise-induced synchronized states. Note that phase locking occurs only when the 
oscillator frequency is approximately resonant to the impulses, whereas the noise-induced synchronization occurs 
almost everywhere. 

Dependence of the Lyapunov exponent (A) on the shape parameter a for several values of oj is plotted in Fig.[5Ja). 
When the natural frequency is resonant to the mean inter- impulse interval (ui = SI = 1/(t) = 1), (A) smoothly 
decreases as we move from the noise- induced synchronization (a = l) to the phase- locking (a = 1000), indicating 
that the phase locking is more stable than the noise-induced synchronization. However, for oscillator frequencies 
near the edges of the phase-locking region (w = 0.91 and 1.09), (A) differs only little between a = 1 and a = 1000, 
so that the noise-induced synchronization can be nearly stable as the phase locking. Thus, though the Lyapunov 
exponent (A) reflects the synchronization dynamics, difference between the two types of synchronization cannot be 
fully characterized just by looking at (A). 

Figure [7|b) shows dependence of the variance of the Lyapunov exponent exponent var(A) on the shape parameter 
a and on the oscillator frequency u, and Figure [5Jb) shows the dependence on the shape parameter a for u — 0.91, 1, 
and 1.09. In contrast to the Lyapunov exponent itself, the variance var(A) always decreases considerably as a is 
increased from the noise-induced synchronization (a = 1) to the phase-locking (a = 1000), because the phase PDF 
becomes much broader for smaller a as we have seen in Fig. [3] Thus, the difference between the phase locking and 
the noise- induced synchronization can be captured more clearly by the fluctuations in the Lyapunov exponent, rather 
than by its long-time average. 

The difference can also be clearly visualized by plotting the ratio of the Lyapunov exponent and its variance 
(A)/var(A) as in Fig. |7fc) and Fig. [SJc). Actually, this ratio characterizes the intermittent dynamics of the phase 
differences between the oscillators as we explain below. 
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C. Intermittent dynamics of phase differences 

Fluctuations of the Lyapunov exponent manifests itself in the dynamics of the phase difference. In Fig. [9j typical 
time sequences of the phase difference A„ between the two oscillators in the synchronized states are shown. Three 
different shape parameters (a = 1, 30, 1000) yielding nearly the same Lyapunov exponent (A) are chosen (with uj = 0.91 
fixed). In noise-induced synchronization (Fig.|9|a)) with small a, the phase difference often exhibits big bursts, because 
the growth rate fluctuates strongly (large var(A)). As the parameter a increases, variation of the phase difference 
becomes smaller (Fig. |U(b)) and, in the phase-locked state, the phase difference stays small and rarely exhibits large 
bursts (Fig.JSfc)), because fluctuations of A are very weak (small var(A)). The peculiar dynamics shown in Fig.[SJa) 
is a direct consequence of the fluctuations of the linear growth rates and intimately related to the well-known on-off 
(or modulational) intermittency [l(| [H, H3] ■ 

Figure [10] shows the stationary PDFs of phase differences between the oscillators for differing values of the shape 
parameter a in normal scales (a) and in log-log scales (b). Reflecting the degree of the fluctuations of A, the width 
of the PDFs varies widely. The tail part of the distribution exhibits power-law decay, which is broader for smaller 
a, indicatin g th at the phase differences exhibit big bursts (phase differences) much more frequently. As shown in 
Ref. 0, Hg-Ull, the exponent of the power-law tail is approximately given by the ratio (A)/var(A) that we plotted 
Fig- Etc)- Thus, degree of intermittency in the synchronization dynamics also characterizes the difference in noise- 
induced synchronization and phase locking clearly. 
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FIG. 9: (Color online) Dynamics of the phase difference A„ at three different values of the shape parameter a. (a) Noise-induced 
synchronization (a = 1). (b) Synchronization induced by intermediate impulses (a — 30). (c) Phase locking (a — 1000). 



VI. INFORMATION ANALYSIS 

We have observed that the joint PDFs of the oscillator phases and the impulse intervals exhibit distinct characteris- 
tics depending on the shape parameter. In this section, to quantify correlations between the oscillator phases and the 
impulse intervals for different types of synchronization dynamics, we introduce information measures. In particular, 
we focus on the mutual dependence or "causality" of a pair of phase oscillators driven by a common impulse source. 
This leads us to physically distinct interpretations of the two synchronization scenarios. 
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FIG. 10: (Color online) Stationary probability distribution P(9 A — 8 B ) of the phase difference 9 A — 8 B of two oscillator calculated 
by direct numerical simulations of the random phase map. Black solid line, blue dashed line and red dotted line plot the results 
for a = 1,30 and 1000, respectively. (a)P(fl' 4 - 6 B ). (h)P{6 A - 6 B ) using logarithmic scaling on the vertical axis. 



A. Information measures 



1. Mutual information 

The mutual information J(X; Y) [39[ is a nonnegative measure that quantifies mutual dependence of two random 
variables A and Y (39}- It is defined as 

J{X;Y) = Y, P(*,v)to-T%fr, (20) 

where p{x, y) is a joint probability distribution of X and Y, and p{x) and p(jj) are probability distributions of X and 
Y, respectively. We often normalize J{X, Y) by its upper limit as 

^ y > = .,.„{ W(y)) - < 21 » 

where H(X) and H(Y) are entropies of X and Y, respectively. The normalized mutual information satisfies < 
j(X; Y) < 1. Assuming, for example, that X represents the inter-impulse interval r and Y the oscillator phase 9 just 
after this interval, J(X; Y) or j(X; Y) quantifies how much the oscillator phase tells us about the impulses. 



2. Interaction information 

The interaction information I(X;Y;Z) [13, EH is a generalization of the mutual information J(X;Y) to three 
random variables A, Y , and Z, and quantifies mutual dependence among them. It is defined as 

I(X;YZ)= > p{x, y In — — — - > p(x, y, z) In — — -, (22 

' p [x )p(v) p(x,z)p(y,z) 

where p(x, y, z) is a joint probability distribution function of A, Y , and Z . I(X; Y; Z) is symmetric with respect to 
permutations of A, Y, and Z. The first term is simply the mutual information J(X;Y) between A and Y, and the 
second term on the right hand side gives conditional mutual information J(A; Y\Z) between A and Y given Z . Thus, 
I(X;Y; Z) can be expressed as 

I(X;Y;Z) = J(X;Y)-J(X;Y\Z), (23) 

which provides us with some intuitive meaning of /(A; Y; Z). Namely, it measures the effect of knowing Z in guessing 
the mutual dependence of A and Y. 

Unlike J(A; Y), I(X; Y; Z) can take both positive and negative values. The sign of /(A; Y; Z) tells us about how 
the three random variables A, Y, and Z depend on each other, as illustrated in Fig. [IT] 
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(i) If I{X\ Y; Z) = 0, at least one of the three random variables is independent of the others. For example, knowing 
Z has no effect in guessing the mutual information between X and Y, J(X;Y\Z) = J(X;Y). All variables are 
mutually dependent if I(X; Y; Z) ^ 0. 

(ii) If I(X; Y; Z) > 0, one of the three variables tends to determine the other two variables. For example, knowing 
Z lowers the information on X gained by observing Y (or vice versa), i.e., J(X;Y\Z) < J(X;Y). In particular, 
the maximum value of I(X;Y; Z), given by mm{H(X),H(Y),H(Z)}, is attained when one of the three variables 
dominates the other two, for example, as X = hx(Z) and Y = hy{Z) where hx and hy are some maps. 

(hi) if I(X;Y;Z) < 0, knowing e.g. Z helps in gaining information on X by observing Y, so that J(X;Y\Z) > 
J(X; Y). In particular, I(X; Y; Z) takes its maximum value when two of the random variables dominate the remaining 
one, for example, when Z = h(X,Y) with some map h and moreover J(Y;Z) = (Y and Z are independent). 
The maximum value is equal to the negative conditional entropy —H(Z\X) of the dominated variable Z given the 
dominating variable X. 

We normalize the interaction information I(X; Y; Z) by its maximum as 

so that < i(X; Y; Z) < 1 is satisfied. In addition to the above normalization, I(X; Y; Z) may also be normalized by 
the mutual information J(X;Y) as 

1 Y > Z) = J{X:Y) = 1 - J(X;Y) ' (25) 

The quantity i'(X;Y; Z) satisfies < i'(X;Y; Z) < 1 (because I(X;Y;Z) > in the present case as shown below), 
quantifies how much information Z conveys about the relationship between X and Y. We use these normalized 
i(X;Y; Z) and i'(X;Y; Z) to quantify mutual dependence among r, 9±, and 02- 
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/(Z;7;Z) = I(X;Y;Z)>0 I(X;Y;Z)<0 

FIG. 11: (Color online) Dependence diagram implied by the interaction information I(X;Y; Z). Each arrow indicates the 
dependence between the variables. 



B. Numerical results 



The mutual information and the interaction information arc calculated from the PDFs obtained by numerically 
solving the Frobenius-Perron equations (fTU)) and (fT2)l . We discretize the inter-impulse intervals and the two oscillator 
phases using 100 grid points and use the resulting coarse-grained discrete probability distributions in the calculation. 
Note that we need the joint probability distribution p(9 A ,9 B 7 r) to calculate the interaction information, which 
corresponds to 100 3 grid points. Such a calculation is only feasible by the Frobenius-Perron approach. 



1. Mutual information J (t; 8) 

Figure I12f a) shows the mutual information J(r; 9) of the inter-impulse interval r and the phase 9 just after this 
interval, as well as its upper limit H(t) (entropy of the inter-impulse intervals) as functions of the shape parameter 
a for the resonant situation, T = l/u> = (r) = 1. H(t) takes its maximum value in the Poisson limit (a = 1), 
monotonously decreases as the shape parameter a is increased, and almost vanishes for nearly periodic impulses 
(a = 10 5 ). The raw mutual information J(t;9) takes a small value for Poisson impulses (a = 1) and gradually 
increases with the shape parameter a. J(r; 9) takes its maximum value at a = am — 1000, then decreases again, and 
almost vanishes when the impulses become nearly periodic (a = 10 5 ). 
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Figure [T2lb) shows the normalized mutual information j(r;9) = J(t;9)/H(t). It increases smoothly from a small 
value (nearly 0) to 1 as the shape parameter a is increased from a = 1 (Poisson) to a = 10 5 (nearly periodic). This 
indicates that the oscillator phase has only little dependence on the inter-impulse interval just before it is measured 
in the Poisson case, whereas the oscillator phase possesses almost complete information about the interval in the 
periodic case. These results are consistent with the shape of the single-oscillator phase PDF shown in Fig. E| 

It is interesting to note that the raw mutual information J(r; 9) has a peak at the intermediate shape parameter a. 
If we regard our model as an information channel, the output oscillator phase conveys information about the input 
impulse interval most efficiently at this value. This can be interpreted as follows. For periodic impulses, even though 
9 is locked to t, the interval r has no information because H{t) = and therefore J(r; 9) = 0. For Poisson impulses, 
though H(t) is relatively large, 9 does not faithfully represent r. As a compromise, J(r; 9) is maximized at the 
intermediate value of a. 




a a 



FIG. 12: (Color online) (a) The mutual information J(8;t) (black solid line) and H(t) (red dashed line), (b) Normalized 
mutual information j(6;r) = J(9;t)/H(t). 



2. Interaction information I(r;6 A ;6 B ) 

More interesting insight can be attained by examining the interaction information I(t;9 a ;9 b ) among the inter- 
impulse interval and the phases of the two driven oscillators. Figure I13f a) shows the interaction information 
I(t;9 a ;9 b ), its upper limit H(t), and mutual informations J(9 A ;9 B ) and J(t;9 a ) as functions of the shape pa- 
rameter a for u} = 1 (resonant). The intermediate peak of the raw interaction information arises due to the same 
reason as in the previous case of the raw mutual information. Interaction informations normalized in three ways, 
h{t-9 a -9 b ) = I(t:9 a ;9 b )/H(t), i 2 { T ;6 A ;0 B ) = I(t;9 a ;9 b )/J(9 a ;9 b ), and i 3 {t;9 a ;9 b ) = I(t;9 a ;9 b )/J(t;9 a ), 
are plotted against a in Fig. [T3l b). 

The normalized interaction information ii(r; 9 A ; 9 B ) decreases as the shape parameter is decreased from a — 50000 
(nearly periodic) to a = 1 (Poisson). For nearly periodic impulses, ii(r;9 A ;9 B ) is nearly 1, implying that r, 9 A , and 
9 B arc significantly dependent on each other in the phase-locking regime. As a is decreased, ii(r;9 A ;9 B ) gradually 
decreases and, for Poisson impulses, ii(r; 9 A ; 9 B ) almost vanishes, indicating that at least one of the 9 A , 9 B , and t is 
almost independent of the others. 

The normalized interaction information ii(r\9 ,9 B ) also decreases as the shape parameter is decreased from a — 
50000 (nearly periodic) to a = 1 (Poisson). For nearly periodic impulses (a = 50000), i2{r;9 A ;9 B ) is nearly 1, which 
implies that the inter-impulse interval r dominates the oscillator phases 9 A and 9 B . As a is decreased, 12(t;9 a ;9 b ) 
decreases gradually, and for Poisson impulses, it becomes very close to 0, namely, r is almost independent of 9 A and 
9 B and thus r tells very little about 9 A and 9 B . 

The last normalized interaction information iy,{T]9 A ;9 B ) also decreases as the shape parameter is decreased from 
a = 50000 (nearly periodic) to a = 1 (Poisson). For nearly periodic impulses (a = 50000), «3(r; 9 A ; 9 B ) is nearly 1, 
which indicates that the phase 9 B has almost complete information about r and 9 A . It means that all t, 9 a and 9 B 
depend on each other strongly. Unlike 12(t;9 a ;9 b ), i3(r;9 A ;9 B ) takes relatively large values around 0.4 even in the 
Poisson limit, which implies that we can get nearly 50% of the information about r and 9 A just by looking at 9 B . 
This is because the oscillators are synchronized even if they are driven by Poisson impulses, and we can learn much 
about 9 A by knowing 9 B . 
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Though the interaction information itself is symmetric with respect to its three arguments, we can obtain further 
insight by taking into account the drive-response configuration of our model, namely, the fact that the inter-impulse 
interval r has a distinct meaning from the other two phase variables 9 A and 9 B . Therefore, for the phase- locked 
situation with ii(r;9 A ;9 B ) ~ 1, we may consider that the phases 9 A and 9 B are predominantly controlled by the 
driving impulse r, as implied from the condition giving the maximum value of 7(r; 9 A ; 9 B ). On the other hand, in the 
Poisson case with ii(r; 9 A ;9 B ) close to 0, we may conclude that t is nearly independent of 9 A and 9 B . This indicates 
that the phases have very little information about the driving impulse in the noise-induced synchronized state. 

The above results lead us to two distinct interpretations of the two synchronization scenarios, as schematically 
summarized in Fig. Qj] in the phase locking, the two oscillators are simply oscillators dominated by the impulses, 
whereas they are almost free from the impulses but still behave synchronously in the noise-induced synchronization. 
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FIG. 13: (Color online) (a) Interaction information I(6 a ;6 b ;t) (black solid line) and normalization factors H(r) (red dashed 
curve), J{6 A ;6 B ) (blue dot-dashed curve) and J(t;6 a ) ( green dotted curve) plotted against the shape parameter a. (b) Nor- 
malized interaction information. i 1 (9 A ; 6 B - r) = I \9 A ; 6 B ; r) / H (r) (red dashed curve), i 2 (6 A ;8 B -r) = I(8 A ;6 B : ;r) / J{8 A ;6 B ) 
(blue dot-dashed curve), and iz(6 A ;6 B ;t) = I(6 A \0 B ;t) / J(t\6 a ) (green dotted curve) 
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FIG. 14: Schematic diagrams showing mutual dependence among the impulse interval and the oscillator phases inferred from 
the interaction information, (a) Noise-induced synchronization, (b) Phase locking. 



VII. SUMMARY 



In the present study, we analyzed uncoupled oscillators driven by gamma impulses that smoothly interpolate between 
periodic and Poisson inter-impulse intervals to quantify the difference in synchronization dynamics between these two 
limiting cases. By examining the dependence of various quantities on the shape parameter of the gamma impulses, 
we have revealed a clear difference between the two types of synchronization. Using phase distributions, frequency 
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detuning, statistics of Lyapunov exponents, and information-theoretic measures, we have quantitatively confirmed our 
original intuition, namely, that the oscillators are principally dominated by the impulses and mutually synchronized 
as its consequence for the periodic driving, whereas the oscillators are not entrained by the impulses (at least not 
strongly influenced by the inter-impulse interval just before the phase is measured) but they still attain coherence in 
the case of Poisson driving. 

Generalization of our results for uncoupled oscillators to mutually interacting oscillators, in particular with delay, 
under common or correlated gamma impulses is an interesting future subject, where not only synchronization but 



qualitatively different behaviors 42j can be expected. More detailed quantification of the input-output and the output- 
output relations of the impulse-driven oscillators from the viewpoint of information transfer would also be interesting, 
where effect of the functional shape of the phase response curve and its optimization would be important issues. As 
we found in the measurement of mutual and interaction information, the transition between the phase locking and 
the noise-induced synchronization may not be simply monotonic from such a viewpoint, implying the importance of 
intermediate randomness of the input impulses. We plan to push forward in these directions in our future studies. 

We thank H. Suctani and H. Hata (Kagoshima), and Y. Tsubo and J. Teramae (RIKEN) for illuminating comments. 
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Appendix. Fourier representation of the Frobenius-Perron equation 



Here we introduce Fourier representation of the Frobenius-Perron equations (|10[) and (|12|). Based on this represen- 
tation, we can use the Fast Fourier Transform (FFT) algorithm to numerically calculate convolution integrals included 
in the equations and drastically reduce calculation costs. 

Firstly, to numerically convolute P n and S, we transform the Dirac delta function in the Frobenius-Perron equation, 
in such a way that it becomes an explicit function of the variable of integration 9' rather than the original phase 
variable 9. 



°° />1 />00 f-C 

Pn+l(0) = V / d& \ 0lT 



dr]W(T)P n (6')R(ri)d (9 - F{9' + 77) -wr-k) (26) 



00 pi pOO pOO 

= V / dff / dr d V W(T)P n (e')R(ri) (27) 



dF-^iO - lot - k) , , 1, . \ 

—6 (6' - F- l {6 -lot -k) + ?/) , (28) 



d6 

where we have defined an inverse function of the phase map, 

F(x) — F~ 1 (x) mod 1. (29) 

Since F is a phase map, 

F~ 1 (9 -ujT-k) =F- 1 {6 - lot) - k = F{9 - lot) + d - k (30) 
holds for F^ 1 (9 — lot) € [d,d+ 1), where d is an integer. Using this, we obtain 



P n+1 (9) = V f dff r dTW{T)P n {6') d[H6 " r) +d k] R (F(9 -uT)+d-k-9') 
= J' d0' J™ dTW{T)P n {6') dHe de UT) f] R(F(9-ojT)-k-9') 



where 



P n+1 {9)= T dTW^T) ^ 9 -^ [ d&Y^v^e^'Y^s^^-^) 
Jo d6 J ^ ^ 



(31) 
(32) 



A:— — 00 

1 d9' P dTW(T)P n (9')^—^S (F(9 - cot) - 9') , (33) 
,/n d9 



S{F{9-lot)-9') = R(T{6-urr)-k-er). (34) 

A;— — 00 

Now, using the discrete Fourier transformation, we can convolute P n and S as 

dF{9-LOT) 



(35) 



' dTW(T)^-^l Y.P^ s ^ m(e ^ T) ( 36 ) 



pOO 

dTW(T)D m {9-L0T), (37) 

Jo 
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where p n ,m and s m are mth Fourier coefBcicnts of P n and 5, respectively: 

p n ,m= [ d6'P n (9')e- 2 ™ 9 ' (38) 



i 

d9'S(6')e- 27rm9 ' , (39) 



o 



and 



D m (6 - cut) = ^"^ e^-"). (40) 
at* 

Next, we need to convolute the functions W(t; a, 6) and D m (9 — lot). From the definition of the gamma distribution, 
Eq. ©, 

W(t; a, b)dr = W{lot; a, Lob)d{uT). (41) 

holds. Using this, we obtain 

/>OC 

Pn+i(0) =y2,p n , m s m / dHW(wr;a,wt)D m (fl-wr) (42) 
Jo 

(43) 



Pn, m s m / d(LUT)2_^Wie 2_^d rn , a e K ' 

m ; Q 

En s IU p 27 ™ " ('44 1 ) 

where u> Q and (i miQ are ath Fourier coefficients of W (lot; a, cob) and D m , respectively: 

/>0O 

W a =/ dW(i;a,w6)e~ 2 ™ f (45) 
Jo 

rf m>a = / d6D m {9)e- 2mae . (46) 
Jo 

Therefore, we can write the evolution of Fourier coefficient of P n (9) as 

Pn-\-\,a ^ Pn,m^mdm,a^^oc ; (^^) 
m 

which can easily be calculated numerically. 

Similarly, we can write the evolution of p n ,m,i, which is the (m, /)th Fourier coefficient of the joint probability 
density function P n (9 A , s ), as 

di,pW a+p , (48) 

ill,! 

which is also easy to calculate numerically. 



